PTEN R Notebook

Multiplex assessment of protein variant abundance by massively parallel sequencing VAMP-seq - multiplex assay that uses fluorescent reporters to measure the steady-state abundance of protein variants in cultured human cells (each cell expresses a single variant directly fused to EGFP…the stability of the variant dictates the abundance of the EGFP fusion and, accordingly, the green fluorescence signal of the cell) - used to assess PTEN and TPMT variants

# Scatter plot of VAMP-seq scores relative to variant position in protein, colored by secondary struct
#essentially getting rid of the last row (wt)
pten1_proc_wt <- pten1_proc[!is.na(pten1_proc$position),]
#creating new column that specifices secondary structure (uses info from helix and sheet columns)
pten1_proc_wt$secondary_struct <- ifelse(is.na(pten1_proc_wt$helix), "unknown",
                        ifelse(pten1_proc_wt$helix==1, "helix",
                        ifelse(pten1_proc_wt$sheet==1, "sheet",
                        ifelse(pten1_proc_wt$helix==0, "neither",
                        "unknown"))))
pten_pos <- ggplot(pten1_proc_wt, aes(x=position, y=score, colour=secondary_struct))+ geom_point(size=.3) + scale_x_continuous(minor_breaks = seq(0, 420, 5)) + scale_color_manual(values=c("#FF4848", "#00C853", "#5757FF", "#A9A9A9")) +ylab("VAMP-seq score")+xlab("Position in PTEN")+labs(colour="Secondary Structure")+ggtitle("PTEN scores in relation to protein structure") + geom_vline(xintercept=27, color="black", size=.1) + geom_vline(xintercept=55, color="black", size=.1) + geom_vline(xintercept=70, color="black", size=.1) + geom_vline(xintercept=85, color="black", size=.1) + geom_vline(xintercept=164.5, color="black", size=.1) + geom_vline(xintercept=212, color="black", size=.1) + geom_vline(xintercept=267.5, color="black", size=.1) + geom_vline(xintercept=343.5, color="black", size=.1)+theme_bw()
#colored by change in hydrophobicity instead of secondary struct
pten_hydro <- ggplot(pten1_proc_wt, aes(x=position, y=score, colour=(hydro2-hydro1)))+ geom_point(size=.3, alpha = 0.3) + scale_x_continuous(minor_breaks = seq(0, 420, 5)) + ylab("VAMP-seq score")+xlab("Position in PTEN")+labs(colour="Hydrophobicity")+ggtitle("PTEN scores in relation to change in hydrophobicity") + geom_vline(xintercept=27, color="black", size=.1) + geom_vline(xintercept=55, color="black", size=.1) + geom_vline(xintercept=70, color="black", size=.1) + geom_vline(xintercept=85, color="black", size=.1) + geom_vline(xintercept=164.5, color="black", size=.1) + geom_vline(xintercept=212, color="black", size=.1) + geom_vline(xintercept=267.5, color="black", size=.1) + geom_vline(xintercept=343.5, color="black", size=.1)
#repeat for tpmt
tpmt1_proc_wt <- tpmt1_proc[!is.na(tpmt1_proc$position),]
tpmt1_proc_wt$secondary_struct <- ifelse(is.na(tpmt1_proc_wt$helix), "unknown",
                        ifelse(tpmt1_proc_wt$helix==1, "helix",
                        ifelse(tpmt1_proc_wt$sheet==1, "sheet",
                        ifelse(tpmt1_proc_wt$helix==0, "neither",
                        "unknown"))))
tpmt_pos <- ggplot(tpmt1_proc_wt, aes(x=position, y=score, colour=secondary_struct))+ geom_point(size=.3) + scale_x_continuous(minor_breaks = seq(0, 405, 5)) + scale_color_manual(values=c("#FF4848", "#00C853", "#5757FF", "#A9A9A9")) +ylab("VAMP-seq score")+xlab("Position in TPMT")+labs(colour="Secondary Structure")+ggtitle("TPMT scores in relation to protein structure") + geom_vline(xintercept=47, color="black", size=.1) + geom_vline(xintercept=78, color="black", size=.1) + geom_vline(xintercept=122.5, color="black", size=.1) + geom_vline(xintercept=140, color="black", size=.1) + geom_vline(xintercept=165, color="black", size=.1) + geom_vline(xintercept=194, color="black", size=.1) + geom_vline(xintercept=209, color="black", size=.1)+theme_bw()
#creating a copy of tpmt1_proc_wt, addition of column that notes change from one secondary struct to another (for violin plot purposes) (plot tpmt_pos_vp to see)
tpmt_colors <- tpmt1_proc_wt
tpmt_colors$fact <- rep(10, nrow(tpmt_colors))
temp <- 1
for(i in 1:(length(tpmt_colors$fact)-1)) {
  if (tpmt_colors$secondary_struct[i] != tpmt_colors$secondary_struct[i+1]) {
    tpmt_colors$fact[i] <- temp
    temp <- temp + 1
  } else {
  tpmt_colors$fact[i] <- temp
  }
}
tpmt_colors$fact[length(tpmt_colors$fact)] <- temp
tpmt_pos_vp <- ggplot(tpmt_colors, aes(x=position, y=score))+ geom_violin(data=tpmt_colors[c(1:2783, 2798:4000),], aes(fill=as.character(fact), colour = factor(TRUE)), draw_quantiles = c(0.5), scale = "width") + 
scale_fill_manual(values=c("1" = "#A9A9A9", "2" = "#00C853", "3" = "#FF4848", "4" = "#00C853","5" = "#FF4848", "6" = "#00C853","7" = "#5757FF", "8" = "#00C853","9" = "#FF4848","10" = "#00C853","11" = "#5757FF", "12" = "#00C853","13" = "#FF4848", "14" = "#00C853", "15" = "#5757FF", "16" = "#00C853", "17" = "#5757FF", "18" = "#00C853", "19" = "#5757FF", "20" = "#00C853", "21" = "#5757FF", "22" = "#00C853", "23" = "#FF4848", "24" = "#5757FF", "25" = "#00C853", "26" = "#FF4848", "27" = "#00C853", "28" = "#5757FF", "29" = "#00C853", "30" = "#5757FF", "31" = "#00C853")) + scale_colour_manual(values = c("black")) + scale_x_continuous(minor_breaks = seq(0, 405, 5)) + ylab("VAMP-seq score")+xlab("Position in TPMT")+labs(colour="Secondary Structure")+ggtitle("TPMT scores in relation to protein structure") + geom_vline(xintercept=47, color="black", size=.1) + geom_vline(xintercept=78, color="black", size=.1) + geom_vline(xintercept=122.5, color="black", size=.1) + geom_vline(xintercept=140, color="black", size=.1) + geom_vline(xintercept=165, color="black", size=.1) + geom_vline(xintercept=194, color="black", size=.1) + geom_vline(xintercept=209, color="black", size=.1)
# Scatter plot with change of hydrophobicity along x-axis, abundance score along y
pten_hydro1 <- ggplot(pten1_proc_wt, aes(y=score, x=(hydro2-hydro1)))+ geom_point(size=0.5, alpha = 0.3) + ylab("Hydrophobicity")+xlab("VAMP-seq score")+ggtitle("PTEN scores in relation to change in hydrophobicity")
# Violin and scatter plots of 
pten_aa_spread <- ggplot(pten1_proc_wt, aes(y=score, x=end)) + geom_violin(draw_quantiles=c(0.5))
pten_aa_spread1 <- ggplot(pten1_proc_wt, aes(y=score, x=end)) + geom_point(size = 0.5, alpha = 0.3, position = position_jitter(width = 0.1, height = NULL))
plot(pten_pos)

plot(tpmt_pos)

#plot(tpmt_pos_vp)
#plot(pten_hydro)
#plot(pten_hydro1)
plot(pten_aa_spread)

plot(pten_aa_spread1)

# Looking for possible correlation between abundance and hydrogen bonding
# preprocessing data (probably unecessary, removes rows with NA in hbond_sum, adds column that once again denotes secondary structure)
pten1_hbond <- pten1_proc[!is.na(pten1_proc$hbond_sum),]
pten1_hbond$secondary_struct <- ifelse(is.na(pten1_hbond$helix), "unknown",
                        ifelse(pten1_hbond$helix==1, "helix",
                        ifelse(pten1_hbond$sheet==1, "sheet",
                        ifelse(pten1_hbond$helix==0, "neither",
                        "unknown"))))
# Scatter of abundance score vs hbond_sum (colored by secondary struct)
pten_plot_hbond <- ggplot(pten1_hbond, aes(x=hbond_sum, y=score, colour=secondary_struct))+ geom_point(alpha=0.4) + ylab("VAMP-seq score")+xlab("DSSP Sum of hydrogen bonds")+ggtitle("PTEN scores in relation to hydrogen bonding") + scale_color_manual(values=c("#FF4848", "#696969", "#5757FF")) + labs(colour="Secondary Structure") + theme_bw()
plot(pten_plot_hbond)

# (uncolored)
pten_plot_hbond1 <- ggplot(pten1_hbond, aes(x=hbond_sum, y=score))+ geom_point(alpha = 0.2) + ylab("VAMP-seq score")+xlab("DSSP Sum of hydrogen bonds")+ggtitle("PTEN scores in relation to hydrogen bonding") + theme_bw()
plot(pten_plot_hbond1)

#repeat for tpmt
tpmt1_hbond <- tpmt1_proc[!is.na(tpmt1_proc$hbond_sum),]
tpmt1_hbond$secondary_struct <- ifelse(is.na(tpmt1_hbond$helix), "unknown",
                        ifelse(tpmt1_hbond$helix==1, "helix",
                        ifelse(tpmt1_hbond$sheet==1, "sheet",
                        ifelse(tpmt1_hbond$helix==0, "neither",
                        "unknown"))))
tpmt_plot_hbond <- ggplot(tpmt1_hbond, aes(x=hbond_sum, y=score, colour=secondary_struct))+ geom_point(alpha=0.4) + ylab("VAMP-seq score")+xlab("DSSP Sum of hydrogen bonds")+ggtitle("PTEN scores in relation to hydrogen bonding") + scale_color_manual(values=c("#FF4848", "#696969", "#5757FF")) + labs(colour="Secondary Structure") + theme_bw()
plot(tpmt_plot_hbond)

tpmt_plot_hbond1 <- ggplot(tpmt1_hbond, aes(x=hbond_sum, y=score))+ geom_point(alpha = 0.2) + ylab("VAMP-seq score")+xlab("DSSP Sum of hydrogen bonds")+ggtitle("TPMT scores in relation to hydrogen bonding") + theme_bw()
plot(tpmt_plot_hbond1)

#less hydrogen bonds ~ higher abundance

#method1 for visualizing tracks (basic, for visualizing in rstudio)
grid.newpage()
grid.draw(rbind(ggplotGrob(tpmt_dssp_schematic), ggplotGrob(tpmt_pos_mean), ggplotGrob(tpmt_heat), size = "last"))

grid.newpage()
grid.draw(rbind(ggplotGrob(pten_dssp_schematic), ggplotGrob(pten_pos_mean), ggplotGrob(pten_heat), size = "last"))

#plotting mean score vs variant changed to 
tpmt_end_sum <- summarySE(tpmt1_proc_wt, measurevar="score", groupvars="end")
tpmt_end_mean <- ggplot(tpmt_end_sum, aes(x=end, y=score)) +
  geom_bar(position=position_dodge(), stat="identity", colour="#999999") +
  geom_errorbar(aes(ymin=score-sd, ymax=score+sd), width=0.001, position=position_dodge()) +
  ylab("mean abundance") + xlab("variant amino acid") + theme_bw()
plot(tpmt_end_mean)

#plotting scores vs 'end' colored by location/secondary structure
tpmt_end_scores_b <- ggplot(data=subset(tpmt1_proc_wt, sheet==1), aes(x=end, y=score, colour=position)) +
  geom_violin(data=subset(tpmt1_proc_wt, sheet==1), draw_quantiles=0.5, scale = "width") +
  geom_point(data=subset(tpmt1_proc_wt, sheet==1), size=.3, alpha = 0.6, position=jitter2) + ylab("abundance") + xlab("variant aa") +
  scale_y_continuous(limits = c(-0.8, 2.2))
tpmt_end_scores_a <- ggplot(data=subset(tpmt1_proc_wt, helix==1), aes(x=end, y=score, colour=position)) +
  geom_violin(data=subset(tpmt1_proc_wt, helix==1), draw_quantiles=0.5, scale = "width") +
  geom_point(data=subset(tpmt1_proc_wt, helix==1), size=.3, alpha = 0.6, position=jitter2) + ylab("abundance") + xlab("variant aa") +
  scale_y_continuous(limits = c(-0.8, 2.2))
tpmt_end_scores_n <- ggplot(data=subset(tpmt1_proc_wt, secondary_struct=="neither"), aes(x=end, y=score, colour=position)) +
  geom_violin(data=subset(tpmt1_proc_wt, secondary_struct=="neither"), draw_quantiles=0.5, scale = "width") +
  geom_point(data=subset(tpmt1_proc_wt, secondary_struct=="neither"), size=.3, alpha = 0.6, position=jitter2) + ylab("abundance") + xlab("variant aa") +
  scale_y_continuous(limits = c(-0.8, 2.2))
plot(tpmt_end_scores_b)

plot(tpmt_end_scores_a)

plot(tpmt_end_scores_n)

#tpmt_end_scores_60 <- ggplot(data=subset(tpmt1_proc_wt, position<=60), aes(x=end, y=score, colour=position)) +
#  geom_violin(data=subset(tpmt1_proc_wt, position<=60), draw_quantiles=0.5, scale = "width") +
#  geom_point(data=subset(tpmt1_proc_wt, position<=60), size=.3, alpha = 0.6, position=jitter2) + ylab("abundance") + xlab("variant aa") +
#  scale_y_continuous(limits = c(-0.8, 2.2))
#tpmt_end_scores_120 <- ggplot(data=subset(tpmt1_proc_wt, position>60 & position<=120), aes(x=end, y=score, colour=position)) +
#  geom_violin(data=subset(tpmt1_proc_wt, position>60 & position<=120), draw_quantiles=0.5, scale = "width") +
#  geom_point(data=subset(tpmt1_proc_wt, position>60 & position<=120), size=.3, alpha = 0.6, position=jitter2) + ylab("abundance") + xlab("variant aa") +
#  scale_y_continuous(limits = c(-0.8, 2.2))
#tpmt_end_scores_180 <- ggplot(data=subset(tpmt1_proc_wt, position>120 & position<=180), aes(x=end, y=score, colour=position)) +
#  geom_violin(data=subset(tpmt1_proc_wt, position>120 & position<=180), draw_quantiles=0.5, scale = "width") +
#  geom_point(data=subset(tpmt1_proc_wt, position>120 & position<=180), size=.3, alpha = 0.6, position=jitter2) + ylab("abundance") + xlab("variant aa") +
#  scale_y_continuous(limits = c(-0.8, 2.2))
#tpmt_end_scores_245 <- ggplot(data=subset(tpmt1_proc_wt, position>180 & position<=245), aes(x=end, y=score, colour=position)) +
#  geom_violin(data=subset(tpmt1_proc_wt, position>180 & position<=245), draw_quantiles=0.5, scale = "width") +
#  geom_point(data=subset(tpmt1_proc_wt, position>180 & position<=245), size=.3, alpha = 0.6, position=jitter2) + ylab("abundance") + xlab("variant aa") +
#  scale_y_continuous(limits = c(-0.8, 2.2))
#plot(tpmt_end_scores_60)
#plot(tpmt_end_scores_120)
#plot(tpmt_end_scores_180)
#plot(tpmt_end_scores_245)

# More of above
pten_a_spread1 <- ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "A"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Alanine variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "A"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_manual(values=twenty_color) + theme_bw()

pten_a_aa <- ggplot(pten1_proc_wt, aes(y=score, x=end, colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "A"), aes(group=end), scale = "width") + xlab("Amino acid") + ggtitle("Alanine variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "A"), aes(x=end, y=score), alpha = 0.75, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

pten_r_spread1 <- ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "R"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Arganine variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "R"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_manual(values=twenty_color) + theme_bw()

pten_r_aa <- ggplot(pten1_proc_wt, aes(y=score, x=end, colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "R"), aes(group=end), scale = "width") + xlab("Amino acid") + ggtitle("Arganine variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "R"), aes(x=end, y=score), alpha = 0.75, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

pten_n_spread1 <- ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "N"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Asparagine variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "N"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_manual(values=twenty_color) + theme_bw()

pten_n_aa <- ggplot(pten1_proc_wt, aes(y=score, x=end, colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "N"), aes(group=end), scale = "width") + xlab("Amino acid") + ggtitle("Asparagine variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "N"), aes(x=end, y=score), alpha = 0.75, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

pten_d_spread1 <- ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "D"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Aspartic Acid variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "D"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_manual(values=twenty_color) + theme_bw()

pten_d_aa <- ggplot(pten1_proc_wt, aes(y=score, x=end, colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "D"), aes(group=end), scale = "width") + xlab("Amino acid") + ggtitle("Aspartic Acid variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "D"), aes(x=end, y=score), alpha = 0.75, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

pten_n_hydrodiff <- ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = (hydro2-hydro1))) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "N"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Asparagine variant abundance scores w/ hydrophobicity change")+ geom_point(data=subset(pten1_proc_wt, start== "N"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_distiller(palette = "Spectral") + theme_bw()

pten_d_hydrodiff <- ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = (hydro2-hydro1))) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "D"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Aspartic Acid variant abundance scores w/ hydrophobicity change")+ geom_point(data=subset(pten1_proc_wt, start== "D"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_distiller(palette = "Spectral") + theme_bw()

plot(pten_a_spread1)
plot(pten_a_aa)
plot(pten_r_spread1)
plot(pten_r_aa)
plot(pten_n_spread1)
plot(pten_n_hydrodiff)
plot(pten_n_aa)
plot(pten_d_spread1)
plot(pten_d_hydrodiff)
plot(pten_d_aa)

ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "P"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Proline variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "P"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_manual(values=twenty_color) + theme_bw()

ggplot(pten1_proc_wt, aes(y=score, x=end, colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "P"), aes(group=end), scale = "width") + xlab("Amino acid") + ggtitle("Proline variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "P"), aes(x=end, y=score), alpha = 0.75, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = (hydro2-hydro1))) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "P"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Proline variant abundance scores w/ hydrophobicity change")+ geom_point(data=subset(pten1_proc_wt, start== "P"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_distiller(palette = "Spectral") + theme_bw()

#graphs for Threonine
ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "T"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("T variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "T"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_manual(values=twenty_color) + theme_bw()

ggplot(pten1_proc_wt, aes(y=score, x=end, colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "T"), aes(group=end), scale = "width") + xlab("Amino acid") + ggtitle("T variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "T"), aes(x=end, y=score), alpha = 0.75, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = (hydro2-hydro1))) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "T"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("T variant abundance scores w/ hydrophobicity change")+ geom_point(data=subset(pten1_proc_wt, start== "T"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_distiller(palette = "Spectral") + theme_bw()

# To plot positions with low stddev (excluding nonsense)

#creating new dataframe which includes mean abundance and variance, and discludes nonsense mutations (nonsense mutations cut proteins short and usually result in very low abundance scores, so I've decided to treat them separately from missense/synonymous mutations)
pten_variance <- summarySE((subset(pten1_data, end != "X")), measurevar="score", groupvars="position", na.rm=TRUE)

# Filtering for positions with enough variants and low variance in variant abundance scores
#adjust N (minimum # variants at a position) and stddev to liking
#5, 0.11
#10, 0.15
pten_variance_filtered <- subset(subset(pten_variance, N > 8 ), sd < 0.11)

#nonsense variant scores are graphed as dots, but are not included in the overlaying violin plots
ggplot(no_nonsense, aes(y=score, x=factor(position), colour = end)) + 
  geom_violin(draw_quantiles=c(0.5), data=subset(no_nonsense, no_nonsense$position %in% pten_variance_filtered$position), aes(group=factor(position)), scale = "width") + 
  xlab("Position in PTEN") + ggtitle("Intolerant and tolerant amino acid positions") +
  geom_point(data=subset(pten1_proc_wt, pten1_proc_wt$position %in% pten_variance_filtered$position), aes(x=factor(position), y=score), alpha = 0.85, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

ggplot(no_nonsense, aes(y=score, x=position, colour = end)) + 
  geom_violin(draw_quantiles=c(0.5), data=subset(no_nonsense, no_nonsense$position %in% pten_variance_filtered$position), aes(group=position%%450), scale = "width") + 
  xlab("Position in PTEN") + ggtitle("Intolerant and tolerant amino acid positions") +
  geom_point(data=subset(pten1_proc_wt, pten1_proc_wt$position %in% pten_variance_filtered$position), aes(x=position, y=score), alpha = 0.85, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

# Repeat of above ,just testing different N and stddev parameter 
pten_variance_filtered1 <- subset(subset(pten_variance, N > 10), sd > 0.25)
ggplot(no_nonsense, aes(y=score, x=factor(position), colour = end)) + 
  geom_violin(draw_quantiles=c(0.5), data=subset(no_nonsense, no_nonsense$position %in% pten_variance_filtered1$position), aes(group=factor(position)), scale = "width") + 
  xlab("Position in PTEN") + ggtitle("Intolerant and tolerant amino acid positions") +
  geom_point(data=subset(pten1_proc_wt, pten1_proc_wt$position %in% pten_variance_filtered1$position), aes(x=factor(position), y=score), alpha = 0.85, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

ggplot(no_nonsense, aes(y=score, x=position, colour = end)) + 
  geom_violin(draw_quantiles=c(0.5), data=subset(no_nonsense, no_nonsense$position %in% pten_variance_filtered1$position), aes(group=position%%450), scale = "width") + 
  xlab("Position in PTEN") + ggtitle("Intolerant and tolerant amino acid positions") +
  geom_point(data=subset(pten1_proc_wt, pten1_proc_wt$position %in% pten_variance_filtered1$position), aes(x=position, y=score), alpha = 0.85, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

chosen <- c(61, 68, 105, 108, 123, 127, 130, 132, 135, 155, 165, 173, 174, 246, 323, 335)
ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + 
  geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, pten1_proc_wt$position %in% chosen), aes(group=factor(position)), scale = "width") + 
  xlab("Position in PTEN") + ggtitle("Clinvar path/likely path variant positions") +
  geom_point(data=subset(pten1_proc_wt, pten1_proc_wt$position %in% chosen), aes(x=factor(position), y=score), alpha = 0.85, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

# New track for heat-mean-variance-secondary-struct plot

#this is how pten_variance was defined: pten_variance <- summarySE((subset(pten1_data, end != "X")), measurevar="score", groupvars="position", na.rm=TRUE)
pten_variance_filtered4 <-subset(pten_variance, N>5)
variance_bar <- ggplot(pten_variance_filtered4, aes(y=sd, x=position)) +
  geom_segment(aes(x=position, xend=position, y=0, yend=sd), color="grey68") +
  geom_point(size=0.5) +
  scale_x_continuous(breaks = seq(0, 403, 20), expand = c(0,0)) + scale_y_continuous(expand = c(0,0))+theme(axis.title.x = element_blank(), axis.text.x = element_blank(), legend.position="none")
plot(variance_bar)

gd=ggplot_gtable(ggplot_build(variance_bar))
maxWidth = grid::unit.pmax(ga$widths, gb$widths, gc$widths, gd$widths)
ga$widths <- as.list(maxWidth)
gb$widths <- as.list(maxWidth)
gc$widths <- as.list(maxWidth)
gd$widths <- as.list(maxWidth)
grid.newpage()

##storing, with specified widths!!
#pdf('pten_tpmt_mean_heat_variance.pdf', width=8, height=6)
#grid.arrange(arrangeGrob(gc,gd,ga,gb,nrow=4,heights=c(.1,.15,.15,.6)))
#dev.off()

# Multiplex paper's Fig2b (PTEN) shows trailing tails for the nonsense and synonymous variant scores
# this is to identify and investigate why variants in this tail are so far from their means
pten1_nonsense <- subset(pten1_proc, class == "nonsense")
tpmt1_nonsense <- subset(tpmt1_proc, class == "nonsense")
pten1_synon <- subset(pten1_proc, class == "synonymous")
tpmt1_synon <- subset(tpmt1_proc, class == "synonymous")
pten1_no_missense <- subset(pten1_proc, class == "synonymous" | class == "nonsense")
ggplot(pten1_nonsense, aes(x=score)) + geom_histogram(binwidth=.01, colour="blue", fill="white") + theme_bw()

#+ geom_density()
ggplot(pten1_synon, aes(x=score)) + geom_histogram(binwidth=.01, colour="red", fill="white") + theme_bw()

ggplot(pten1_proc_wt, aes(x=score)) + geom_histogram(data=subset(pten1_proc_wt,class == "nonsense"), fill = "red", alpha = 0.5, binwidth=.01) + geom_histogram(data=subset(pten1_proc_wt,class == "synonymous"), fill = "blue", alpha = 0.5, binwidth=.01) + geom_histogram(data=subset(pten1_proc_wt,class == "missense"), fill = "green", alpha = 0.2, binwidth=.01) + theme_bw()

ggplot(pten1_no_missense, aes(x=score)) + geom_histogram(data=subset(pten1_no_missense,class == "nonsense"), fill = "red", alpha = 0.5, binwidth=.01) + geom_histogram(data=subset(pten1_no_missense,class == "synonymous"), fill = "blue", alpha = 0.5, binwidth=.01) + theme_bw()

# TPMT
#ggplot(tpmt1_nonsense, aes(x=score)) + geom_histogram(binwidth=.01, colour="blue", fill="white") + theme_bw()
#ggplot(tpmt1_synon, aes(x=score)) + geom_histogram(binwidth=.01, colour="red", fill="white") + theme_bw()
# we can see the paper used a different scale, which was ... deceiving
# not many variants are nonsense or synonymous, so the tail is not actually that large for both of them

# Taking a subset of the data, all variants above a certain score (for nonsense), and below a certain score (for synonymous)
#0.55
nonsense_tail <- subset(pten1_nonsense, score > 0.6)
synon_tail <- subset(pten1_synon, score < 0.6)
nonsense_tail$secondary_struct <- ifelse(is.na(nonsense_tail$helix), "unknown",
                        ifelse(nonsense_tail$helix==1, "helix",
                        ifelse(nonsense_tail$sheet==1, "sheet",
                        ifelse(nonsense_tail$helix==0, "neither",
                        "unknown"))))
synon_tail$secondary_struct <- ifelse(is.na(synon_tail$helix), "unknown",
                        ifelse(synon_tail$helix==1, "helix",
                        ifelse(synon_tail$sheet==1, "sheet",
                        ifelse(synon_tail$helix==0, "neither",
                        "unknown"))))
n_tail <- nonsense_tail[,c(1,2,7,30,127)]
s_tail <- synon_tail[,c(1,2,7,30,127)]
n_tail$bp_pos <- (n_tail$position-1)*3
s_tail$bp_pos <- (s_tail$position-1)*3
n_tail

s_tail

# Graphing the positions of these variants in the protein
#just in case there is a discernible pattern
s_tail_pos <- ggplot(s_tail, aes(x=position, y=score, colour=secondary_struct))+ geom_point(size=.3) + scale_x_continuous(minor_breaks = seq(0, 405, 5)) + scale_color_manual(values=c("#FF4848", "#00C853", "#5757FF", "#A9A9A9")) +ylab("VAMP-seq score")+xlab("Position in PTEN")+labs(colour="Secondary Structure")+ggtitle("PTEN synonymous variant tail scores in relation to protein structure") + theme_bw()
plot(s_tail_pos)

n_tail_pos <- ggplot(n_tail, aes(x=position, y=score, colour=secondary_struct))+ geom_point(size=.3) + scale_x_continuous(minor_breaks = seq(0, 405, 5)) + scale_color_manual(values=c("#FF4848", "#00C853", "#5757FF", "#A9A9A9")) +ylab("VAMP-seq score")+xlab("Position in PTEN")+labs(colour="Secondary Structure")+ggtitle("PTEN nonsense variant tail scores in relation to protein structure") + theme_bw()
plot(n_tail_pos)

TPMT_abun_CADD <- ggplot(tpmt_merge, aes(x=abundance_class, y=CADD_raw_rankscore)) + geom_violin(draw_quantiles = c( 0.5))+ylab("CADD raw rankscore")+xlab("Abundance Class") + theme_bw()
plot(TPMT_abun_CADD)

TPMT_abun_SIFT_conv <- ggplot(tpmt_merge, aes(x=abundance_class, y=as.numeric(SIFT_converted_rankscore))) + geom_violin(draw_quantiles = c(0.5))+ylab("SIFT conv rankscore")+xlab("Abundance Class") + theme_bw()
plot(TPMT_abun_SIFT_conv)

TPMT_abun_POLY <- ggplot(tpmt_merge, aes(x=abundance_class, y=as.numeric(Polyphen2_HDIV_rankscore))) + geom_violin(draw_quantiles = c( 0.5))+ylab("Polyphen2 HDIV rankscore")+xlab("Abundance Class") + theme_bw()
plot(TPMT_abun_POLY)

TPMT_abun_POLY1 <- ggplot(tpmt_merge, aes(x=abundance_class, y=as.numeric(Polyphen2_HVAR_rankscore))) + geom_violin(draw_quantiles = c( 0.5))+ylab("Polyphen2 HVAR rankscore")+xlab("Abundance Class") + theme_bw()
plot(TPMT_abun_POLY1)

# Bar plots of what Polyphen and SIFT scores make up each abundance class!
Pred_abun_SIFT <- ggplot(tpmt_merge, aes(abundance_class)) + geom_bar(aes(fill = SIFT_pred)) + ggtitle("Abundance class vs SIFT prediction of Damaging or Tolerated") + theme_bw()
plot(Pred_abun_SIFT)

trial_sep <- tpmt_merge[c(21,23,24,26)]
tpmt_merge_expand <- separate_rows(tpmt_merge, c("Polyphen2_HDIV_score", "Polyphen2_HDIV_pred", "Polyphen2_HVAR_score", "Polyphen2_HVAR_pred"))
Pred_abun_HVAR <- ggplot(tpmt_merge_expand, aes(abundance_class)) + geom_bar(aes(fill = Polyphen2_HVAR_pred)) + ggtitle("Abundance class vs Polyphen2 HVAR predictions") + labs(subtitle = "D: Probably Damaging, P: Possibly Damaging, B: Benign") + theme_bw()
plot(Pred_abun_HVAR)

# Focusing on individual "end", or variant, amino acids! A lot more can be done... and a lot of what I have here is unfocused and too busy to be helpful
twenty_color1 = c("#D02028", "#A4C33B","#53958B", "#E6A3B4", "#C5A0CA", "#554DA0", "#99247E", "#402059", "#82421B", "#7E807E", 'black', "#EDD941", "#F2F08E", "#EEC898", "#E1A12F", "#76C158",  "#BCDDAE", "#85782E", "#315935", "#A1DAE0", "#486EB6")
pten_dssp_schematic1 <- ggplot() +
  geom_segment(aes(x = 1, y = 0, xend = max(pten_extra$position)), yend = 0, size = 1, color = "grey70") +
  geom_point(data = subset(pten_extra, !is.na(xca)), aes(x = position, y = 0), color = "black", size = 1.8) +
  geom_point(data = subset(pten_extra, sheet == 1), aes(x = position, y = 0), color = "pink", size = 1.5) +
  geom_point(data = subset(pten_extra, helix == 1), aes(x = position, y = 0), color = "cyan", size = 1.5) +
  scale_x_continuous(breaks = seq(0, 403, 20), expand = c(0,0)) +
  scale_y_continuous(breaks = NULL, expand = c(0,0)) +xlab("Position in PTEN") + ylab("\n \n \n") +
  theme(panel.border = element_blank(), axis.text.y = element_blank())
aas <- c("A", "C", "P", "X")
aas1 <- c("S", "C", "P", "X")
aas2 <- c("A", "C", "D", "E", "F", "G", "H", "I", "K", "L", "M", "N", "P", "Q", "R", "S", "T", "V", "W", "Y")
pten_a_pos <- ggplot(data=subset(pten1_proc_wt, end=="A"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="A" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Alanine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
#plot(pten_a_pos)
pten_s_pos <- ggplot(data=subset(pten1_proc_wt, end=="S"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="S" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Serine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
pten_c_pos <- ggplot(data=subset(pten1_proc_wt, end=="C"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="C" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Cysteine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
pten_x_pos <- ggplot(data=subset(pten1_proc_wt, end=="X"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="X" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of nonsense variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_a_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

grid.newpage()
grid.draw(rbind(ggplotGrob(pten_s_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

grid.newpage()
grid.draw(rbind(ggplotGrob(pten_c_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

grid.newpage()
grid.draw(rbind(ggplotGrob(pten_x_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

# More of above
#snake at the end (arches over grey region)
pten_p_pos <- ggplot(data=subset(pten1_proc_wt, end=="P"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="P" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Proline variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_p_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

# More of above
#matching snake at the end!
pten_g_pos <- ggplot(data=subset(pten1_proc_wt, end=="G"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="G" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Glycine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_g_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

pten_h_pos <- ggplot(data=subset(pten1_proc_wt, end=="H"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="H" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Histidine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_h_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

# notice 2nd from last grey region
pten_w_pos <- ggplot(data=subset(pten1_proc_wt, end=="W"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="W" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Typtophan variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_w_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))